# total random sites to create
tot <- nrow(restdat)
id <- stri_rand_strings(tot, 4)
dts <- range(restdat$date)
# rnorm(10 * tot, mean(reststat$lon), sd(reststat$lon))
# rnorm(10 * tot, mean(reststat$lat), sd(reststat$lat))
ext <- bbox(tbpoly)
lon <- runif(10 * tot, ext[1, 1], ext[1, 2])
lat <- runif(10 * tot, ext[2, 1], ext[2, 2])
tmp <- SpatialPoints(cbind(lon, lat),
proj4string = crs(tbpoly)
) %>%
.[tbpoly, ] %>%
.[sample(1:nrow(.@coords), tot, replace = F), ]
restdat_rnd <- tibble(id) %>%
mutate(
date = sample(seq(dts[1], dts[2]), tot, replace = T),
top = sample(c('hab', 'wtr'), tot, replace = T)
)
reststat_rnd <- tibble(id) %>%
mutate(
lat = tmp$lat,
lon = tmp$lon
)
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'
# base map
ext <- make_bbox(reststat_rnd$lon, reststat_rnd$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 10, maptype = "toner-lite")
pbase <- ggmap(map) +
theme_bw() +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
# map by restoration type
pbase +
geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21)
# map by date
pbase +
geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = factor(date)), size = 4, pch = 21)
# barplot of date counts
toplo <- restall_rnd %>%
group_by(date)
ggplot(restall_rnd, aes(x = factor(date))) +
geom_bar() +
coord_flip() +
theme_bw() +
theme(
axis.title.y = element_blank()
) +
scale_y_discrete(expand = c(0, 0))
wqmtch_rnd <- get_clo(restdat_rnd, reststat_rnd, wqstat, resgrp = 'top', mtch = mtch)
save(wqmtch_rnd, file = 'data/wqmtch_rnd.RData', compress = 'xz')
head(wqmtch_rnd)
## # A tibble: 6 x 5
## stat resgrp rnk id dist
## <int> <chr> <int> <chr> <dbl>
## 1 47 wtr 1 cTUY 3806.092
## 2 47 wtr 2 Ryvp 6696.894
## 3 47 wtr 3 fwaq 7272.233
## 4 47 wtr 4 0lGd 7962.964
## 5 47 wtr 5 peUI 9097.059
## 6 47 wtr 6 Kf3O 9290.502
##
# plots
# combine lat/lon for the plot
toplo <- wqmtch_rnd %>%
left_join(wqstat, by = 'stat') %>%
left_join(reststat_rnd, by = 'id') %>%
rename(
`Restoration\ngroup` = resgrp,
`Distance (dd)` = dist
)
# restoration project grouping column
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'
# extent
ext <- make_bbox(wqstat$lon, wqstat$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 12, maptype = "toner-lite")
# base map
pbase <- ggmap(map) +
theme_bw() +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank()
) +
geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21) +
geom_point(data = wqstat, aes(x = lon, y = lat), size = 2)
# closest
toplo1 <- filter(toplo, rnk %in% 1)
pbase +
geom_segment(data = toplo1, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)
# closest five percent
fvper <- max(toplo$rnk) %>%
`*`(0.2) %>%
ceiling
toplo2 <- filter(toplo, rnk %in% c(1:fvper))
pbase +
geom_segment(data = toplo2, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)
# closest all combo
toplo3 <- toplo
pbase +
geom_segment(data = toplo3, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)
Get weighted average of project type, treatment (before, after) of salinity for all wq station, restoration site combinations.
salchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf, chgout = TRUE)
salchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf)
save(salchgout_rnd, file = 'data/salchgout_rnd.RData')
save(salchg_rnd, file = 'data/salchg_rnd.RData')
head(salchgout_rnd)
## # A tibble: 6 x 3
## # Groups: stat [2]
## stat cmb cval
## <int> <chr> <dbl>
## 1 6 hab_aft 24.31485
## 2 6 hab_bef 23.64315
## 3 6 wtr_aft 24.55549
## 4 6 wtr_bef 23.84958
## 5 7 hab_aft 24.81490
## 6 7 hab_bef 24.42562
head(salchg_rnd)
## # A tibble: 6 x 4
## stat hab wtr cval
## <int> <fctr> <fctr> <dbl>
## 1 6 hab_aft wtr_aft 24.43517
## 2 6 hab_aft wtr_bef 24.08221
## 3 6 hab_bef wtr_aft 24.09932
## 4 6 hab_bef wtr_bef 23.74636
## 5 7 hab_aft wtr_aft 24.86153
## 6 7 hab_aft wtr_bef 24.65012
Get conditional probability distributions for the restoration type, treatment effects, salinity as first child node in network.
wqcdt_rnd <- get_cdt(salchg_rnd, 'hab', 'wtr')
head(wqcdt_rnd)
## # A tibble: 4 x 5
## hab wtr data crv prd
## <fctr> <fctr> <list> <list> <list>
## 1 hab_aft wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 2 hab_aft wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 3 hab_bef wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 4 hab_bef wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
Discretization of salinity conditional probability distributions:
salbrk_rnd <- get_brk(wqcdt_rnd, qts = c(0.33, 0.66), 'hab', 'wtr')
salbrk_rnd
## # A tibble: 8 x 5
## hab wtr qts brk clev
## <fctr> <fctr> <dbl> <dbl> <dbl>
## 1 hab_aft wtr_aft 26.39854 0.4506171 1
## 2 hab_aft wtr_aft 29.77876 0.8903498 2
## 3 hab_aft wtr_bef 26.19278 0.4321444 1
## 4 hab_aft wtr_bef 29.82180 0.8790707 2
## 5 hab_bef wtr_aft 26.35804 0.4415560 1
## 6 hab_bef wtr_aft 29.83229 0.8784237 2
## 7 hab_bef wtr_bef 26.15228 0.4249726 1
## 8 hab_bef wtr_bef 29.87533 0.8678202 2
A plot showing the breaks:
toplo <- dplyr::select(wqcdt_rnd, -data, -crv) %>%
unnest
ggplot(toplo, aes(x = cval, y = cumest)) +
geom_line() +
geom_segment(data = salbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
geom_segment(data = salbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
facet_grid(hab ~ wtr) +
theme_bw()
Get conditional probability distributions for the restoration type, treatment effects, salinity levels, chlorophyll as second child node in network.
# get chlorophyll changes
chlchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf, chgout = TRUE)
chlchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf)
save(chlchgout_rnd, file = 'data/chlchgout_rnd.RData')
save(chlchg_rnd, file = 'data/chlchg_rnd.RData')
# merge with salinity, bet salinity levels
salbrk_rnd <- salbrk_rnd %>%
group_by(hab, wtr) %>%
nest(.key = 'levs')
allchg_rnd <- full_join(chlchg_rnd, salchg_rnd, by = c('hab', 'wtr', 'stat')) %>%
rename(
salev = cval.y,
cval = cval.x
) %>%
group_by(hab, wtr) %>%
nest %>%
left_join(salbrk_rnd, by = c('hab', 'wtr')) %>%
mutate(
sallev = pmap(list(data, levs), function(data, levs){
# browser()
out <- data %>%
mutate(
saval = salev,
salev = cut(salev, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
salev = as.character(salev)
)
return(out)
})
) %>%
dplyr::select(-data, -levs) %>%
unnest
salchg_rnd <- dplyr::select(allchg_rnd, stat, hab, wtr, salev, saval)
save(salchg_rnd, file = 'data/salchg_rnd.RData', compress = 'xz')
chlcdt_rnd <- get_cdt(allchg_rnd, 'hab', 'wtr', 'salev')
save(chlcdt_rnd, file = 'data/chlcdt_rnd.RData', compress = 'xz')
chlbrk_rnd <- get_brk(chlcdt_rnd, c(0.33, 0.66), 'hab', 'wtr', 'salev')
chlbrk_rnd %>%
print(n = nrow(.))
## # A tibble: 24 x 6
## hab wtr salev qts brk clev
## <fctr> <fctr> <chr> <dbl> <dbl> <dbl>
## 1 hab_aft wtr_aft lo 12.385300 0.4647431 1
## 2 hab_aft wtr_aft lo 17.290408 0.9175791 2
## 3 hab_aft wtr_aft md 7.934825 0.4577565 1
## 4 hab_aft wtr_aft md 10.001101 0.8997290 2
## 5 hab_aft wtr_aft hi 4.118044 0.2745101 1
## 6 hab_aft wtr_aft hi 5.023006 0.7254158 2
## 7 hab_aft wtr_bef lo 14.069949 0.5312290 1
## 8 hab_aft wtr_bef lo 19.673450 0.9323183 2
## 9 hab_aft wtr_bef md 8.242412 0.3989296 1
## 10 hab_aft wtr_bef md 10.204171 0.8325231 2
## 11 hab_aft wtr_bef hi 4.205353 0.2998067 1
## 12 hab_aft wtr_bef hi 4.942181 0.7312806 2
## 13 hab_bef wtr_aft lo 11.636849 0.3631110 1
## 14 hab_bef wtr_aft lo 15.729469 0.8169240 2
## 15 hab_bef wtr_aft md 8.187762 0.4362466 1
## 16 hab_bef wtr_aft md 10.192796 0.8625787 2
## 17 hab_bef wtr_aft hi 4.049786 0.2586025 1
## 18 hab_bef wtr_aft hi 4.789796 0.7109169 2
## 19 hab_bef wtr_bef lo 12.461799 0.3684596 1
## 20 hab_bef wtr_bef lo 17.357014 0.8163132 2
## 21 hab_bef wtr_bef md 9.582194 0.4428301 1
## 22 hab_bef wtr_bef md 12.651273 0.8677568 2
## 23 hab_bef wtr_bef hi 4.137095 0.2927909 1
## 24 hab_bef wtr_bef hi 4.708972 0.7124084 2
Final combinations long format:
chlbar_rnd <- chlbrk_rnd %>%
group_by(hab, wtr, salev) %>%
nest %>%
mutate(
data = map(data, function(x){
brk <- x$brk
out <- data.frame(
lo = brk[1], md = brk[2] - brk[1], hi = 1 - brk[2]
)
return(out)
})
) %>%
unnest %>%
gather('chllev', 'chlval', lo:hi) %>%
mutate(
salev = factor(salev, levels = c('lo', 'md', 'hi')),
chllev = factor(chllev, levels = c('lo', 'md', 'hi'))
)
save(chlbar_rnd, file = 'data/chlbar_rnd.RData', compress = 'xz')
chlbar_rnd %>%
print(n = nrow(.))
## # A tibble: 36 x 5
## hab wtr salev chllev chlval
## <fctr> <fctr> <fctr> <fctr> <dbl>
## 1 hab_aft wtr_aft lo lo 0.46474308
## 2 hab_aft wtr_aft md lo 0.45775646
## 3 hab_aft wtr_aft hi lo 0.27451009
## 4 hab_aft wtr_bef lo lo 0.53122896
## 5 hab_aft wtr_bef md lo 0.39892956
## 6 hab_aft wtr_bef hi lo 0.29980671
## 7 hab_bef wtr_aft lo lo 0.36311099
## 8 hab_bef wtr_aft md lo 0.43624662
## 9 hab_bef wtr_aft hi lo 0.25860246
## 10 hab_bef wtr_bef lo lo 0.36845964
## 11 hab_bef wtr_bef md lo 0.44283010
## 12 hab_bef wtr_bef hi lo 0.29279089
## 13 hab_aft wtr_aft lo md 0.45283601
## 14 hab_aft wtr_aft md md 0.44197255
## 15 hab_aft wtr_aft hi md 0.45090572
## 16 hab_aft wtr_bef lo md 0.40108937
## 17 hab_aft wtr_bef md md 0.43359351
## 18 hab_aft wtr_bef hi md 0.43147385
## 19 hab_bef wtr_aft lo md 0.45381302
## 20 hab_bef wtr_aft md md 0.42633209
## 21 hab_bef wtr_aft hi md 0.45231442
## 22 hab_bef wtr_bef lo md 0.44785357
## 23 hab_bef wtr_bef md md 0.42492674
## 24 hab_bef wtr_bef hi md 0.41961753
## 25 hab_aft wtr_aft lo hi 0.08242091
## 26 hab_aft wtr_aft md hi 0.10027099
## 27 hab_aft wtr_aft hi hi 0.27458419
## 28 hab_aft wtr_bef lo hi 0.06768166
## 29 hab_aft wtr_bef md hi 0.16747694
## 30 hab_aft wtr_bef hi hi 0.26871943
## 31 hab_bef wtr_aft lo hi 0.18307600
## 32 hab_bef wtr_aft md hi 0.13742130
## 33 hab_bef wtr_aft hi hi 0.28908311
## 34 hab_bef wtr_bef lo hi 0.18368680
## 35 hab_bef wtr_bef md hi 0.13224315
## 36 hab_bef wtr_bef hi hi 0.28759158
Discretesize chlorophyll data, all stations:
# discretize all chl data by breaks
chlbrk_rnd <- chlbrk_rnd %>%
group_by(hab, wtr, salev) %>%
nest(.key = 'levs')
allchg_rnd <- allchg_rnd %>%
group_by(hab, wtr, salev) %>%
nest %>%
full_join(chlbrk_rnd, by = c('hab', 'wtr', 'salev')) %>%
mutate(
lev = pmap(list(data, levs), function(data, levs){
out <- data %>%
mutate(
lev = cut(cval, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
lev = as.character(lev)
)
return(out)
})
) %>%
dplyr::select(-data, -levs) %>%
unnest %>%
rename(
chlev = lev,
chval = cval
)
save(allchg_rnd, file = 'data/allchg_rnd.RData', compress = 'xz')
A bar plot of splits:
ggplot(chlbar_rnd, aes(x = chllev, y = chlval, group = salev, fill = salev)) +
geom_bar(stat = 'identity', position = 'dodge') +
facet_grid(hab ~ wtr) +
theme_bw()
A plot showing the breaks:
toplo <- dplyr::select(chlcdt_rnd, -data, -crv) %>%
unnest %>%
mutate(
salev = factor(salev, levels = c('lo', 'md', 'hi'))
)
chlbrk_rnd <- unnest(chlbrk_rnd)
ggplot(toplo, aes(x = cval, y = cumest, group = salev, colour = salev)) +
geom_line() +
geom_segment(data = chlbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
geom_segment(data = chlbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
facet_grid(hab ~ wtr, scales = 'free_x') +
theme_bw()